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ABSTRACT 

We present results for Vela C obtained during the 2012 flight of the Balloon-borne Large Aperture 
Submillimeter Telescope for Polarimetry (BLASTPol). We mapped polarized intensity across almost 
the entire extent of this giant molecular cloud, in bands centered at 250, 350, and 500/xm. In this 
initial paper, we show our 500/xm data smoothed to a resolution of 2^5 (approximately 0.5pc). We 
show that the mean level of the fractional polarization p and most of its spatial variations can be 
accounted for using an empirical three-parameter power-law ht, poc where N is the 

hydrogen column density and S is the polarization-angle dispersion on 0.5 pc scales. The decrease of 
p with increasing S is expected because changes in the magnetic field direction within the cloud volume 
sampled by each measurement will lead to cancellation of polarization signals. The decrease of p with 
increasing N might be caused by the same effect, if magnetic field disorder increases for high column 
density sightlines. Alternatively, the intrinsic polarization efficiency of the dust grain population 
might be lower for material along higher density sightlines. We find no significant correlation between 
N and S. Comparison of observed submillimeter polarization maps with synthetic polarization maps 
derived from numerical simulations provides a promising method for testing star formation theories. 
Realistic simulations should allow for the possibility of variable intrinsic polarization efficiency. The 
measured levels of correlation among p, N^ and S provide points of comparison between observations 
and simulations. 

Subject headings: instrumentation: polarimeters, ISM: dust, extinction, ISM: magnetic ffelds, ISM: 
individual objects (Vela C), stars: formation, techniques: polarimetric 
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1. INTRODUCTION 

The Balloon-borne Large Aperture Submillimeter 
Telescope for Polarimetry (BLASTPol; Galitzki et al. 
2014) is sensitive to magnetic field structure ranging from 
scales of entire giant molecular clouds (GMCs) down to 
cores (for nearby clouds). In this paper we present a very 
sensitive survey of the star-forming region Vela C from 
the 2012 flight of BLASTPol. Our goal is to quantify the 
dependence of polarization fraction on column density, 
temperature, and local magnetic field disorder, in order 
to provide empirical formulae that can be used to test 
numerical simulations of molecular clouds. These obser¬ 
vations are timely because the role that magnetic fields 
play in the formation of molecular clouds and their sub¬ 
structures persists as an outstanding question in the un¬ 
derstanding of the detailed mechanics of star formation 
(McKee & Ostriker 2007). Strong magnetic fields that 
are well coupled to the gas can inhibit or slow down grav¬ 
itational collapse of gas in the direction perpendicular 
to the cloud magnetic field lines (Mouschovias & Giolek 
1999). This in turn can contribute to the low observed 
star formation efficiency seen in molecular clouds. Nu¬ 
merical simulations of molecular clouds show that mag¬ 
netized clouds differ from unmagnetized clouds in cloud 
density contrasts (Kowal et al. 2007) and in star forma¬ 
tion efficiency (Myers et al. 2014). However, obtaining 
detailed measurements of magnetic fields in molecular 
clouds over a wide range of relevant spatial and density 
scales remains challenging. 

Zeeman splitting in molecular lines can be used to mea¬ 
sure the component of magnetic field strength parallel to 
the line of sight directly (Grutcher 2012). However this 
technique is challenging as the Doppler broadening of 
molecular lines is generally much larger than the Zeeman 
splitting. After many years of careful observations there 
are now several dozen detections of Zeeman splitting in 
molecular lines that trace dense gas (Grutcher 2012). 

An alternative method for studying magnetic fields in 
molecular clouds is to use polarization maps to infer the 
orientation of the magnetic field projected on the plane 
of the sky (4)). Dust grains are believed to align with 
their long axes on average perpendicular to the local 
magnetic field (see Lazarian 2007 for a review). Gurrent 
evidence suggests that radiative torques (RATs) from 
anisotropic radiation fields might be the dominant align¬ 
ment mechanism (Lazarian & Hoang 2007; Andersson 
et al. 2015). Optical and near-IR light from stars that 
passes through a foreground cloud of aligned dust grains 
becomes polarized parallel to 4>. This method has long 
been used to study the magnetic field orientation in the 
diffuse ISM (Hall 1949; Hiltner 1949; Heiles 2000), but 
is not easily applicable for high extinction cloud sight¬ 
lines. However, dust grains emit radiation preferentially 
polarized parallel to their long axes, so that the resulting 
far-infrared/submillimeter thermal emission is polarized 
orthogonal to 4> (Hildebrand 1988). The emission is gen¬ 
erally optically thin. 

The fraction of dust emission that is polarized (p), 
does not give any direct estimate of the magnetic field 
strength. However, it can encode information about the 
dust grain shape and alignment efficiency, angle of the 
field with respect to the line of sight, and changes in 
field direction. Hildebrand (1988) reviews the factors 


that affect p for optically-thin thermal emission from a 
population of grains. First, consider the case of perfect 
spinning alignment of an ensemble of identical grains in a 
uniform magnetic field oriented orthogonally to the line 
of sight (7 = 0°). In this case p will be determined 
by the grains’ optical constants and shape (e.g., ratio 
of axes for the case of oblate spheroids). Next, if the 
grain spin axes are not all exactly parallel to the am¬ 
bient field, the polarization will be reduced by what is 
known as the Rayleigh reduction factor (Greenberg 1968; 
Lazarian 2007). For this paper we define the “intrinsic 
polarization efficiency” as the polarization p of the emis¬ 
sion from such an ensemble of imperfectly aligned grains. 
The measured polarization fraction can be less than this 
intrinsic polarization efficiency if there are variations in 
magnetic field direction within the conical volume being 
sampled by an observation. Finally, for arbitrary values 
of 7 , the polarization is proportional to cos^( 7 ). 

Comparisons between statistical properties of observed 
polarization maps and synthetic observations of 3-D nu¬ 
merical models of star formation are a promising method 
for constraining molecular cloud physics. Examples in¬ 
clude Falceta-Gongalves et al. (2008) as well as his¬ 
tograms of relative orientations (HRO, Soler et al. 2013; 
Planck Collaboration Int. XXXV 2016). If the intrinsic 
polarization efficiency varies within the cloud, then mea¬ 
surements of the inferred magnetic field orientation will 
be weighted toward the field orientation in regions along 
the line of sight where the intrinsic polarization efficiency 
is high. To use polarization observations to constrain the 
structure of the magnetic field in star-forming clouds it 
is therefore important to understand how the intrinsic 
polarization efficiency varies within molecular clouds. 

The VelaC GMC was discovered by Murphy & May 
(1991) via CO observations of a larger scale structure 
known as the Vela Molecular Ridge. VelaC was later ob¬ 
served in the submillimeter by Netterfield et al. (2009) 
and was found to be a cool molecular cloud in an early 
evolutionary state. At a distance of 700 ± 200 pc (Liseau 
et al. 1992), the cloud subtends 3° on the sky (35pc), 
and contains a large quantity of dense gas (M « 5 x 
10"^ Mq as traced by C^®0 1-0 observations from Ya- 
maguchi et al. 1999). A HerscheP^ survey of Vela C 
by Hill et al. (2011) showed that the cloud could be di¬ 
vided into five subregions at an Ay = 7 mag threshold as 
shown in Figure 1. These subregions show a range of 
cloud substructures, from the apparently cold network 
of overlapping filaments in the South-Nest subregion, to 
the high mass Centre-Ridge, which contains a compact 
HII region, RCW36. 

This paper presents an overview of the BLAST¬ 
Pol 500 ^m maps of Vela C from the 2012 flight. BLAST¬ 
Pol polarization data at 250 ^m and 350 pm are discussed 
in a separate paper on the polarization spectrum of 
Vela C (Gandilo et al. 2016). In Section 2 we describe the 
BLASTPol telescope and BLASTPol 2012 science flight. 
Section 3 gives an overview of the data analysis pipeline 
and Section 4 presents the BLASTPol 500 pm polariza¬ 
tion maps. For comparison with the BLASTPol polariza¬ 
tion data we used spectral energy distribution fits to the 

Herschel is an ESA space observatory with science instru- 
ments provided by European-led Principal Investigator consortia 
and with important participation from NASA. 
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well-calibrated, higher resolution Herschel SPIRE and 
PACS maps to produce maps of Vela C column density 
(V) and dust temperature (T) as described in Section 5. 
We then examine the correlations between the polariza¬ 
tion fraction p and N and T in Section 6, and develop 
a two-variable power-law model of p as a function of N 
and the local polarization-angle dispersion S in Section 
7. Finally, in Section 8, we discuss the implications of 
our power-law model and we place a rough upper limit 
on the degree to which reduced intrinsic polarization ef¬ 
ficiency at high N might bias our polarization measure¬ 
ments toward lower density cloud regions. Our hndings 
are summarized in Section 9. 

2. OBSERVATIONS 

BLASTPol is a high altitude balloon-borne polarime- 
ter that utilizes a 1.8-m diameter aluminium parabolic 
primary mirror, and a 40-cm aluminum secondary mir¬ 
ror. Incoming light is directed onto a series of reimaging 
optics cooled to 1K in a liquid nitrogen-helium cryostat 
(Galitzki et al. 2014). A series of dichroic filters direct 
light onto focal-plane arrays of bolometers (cooled be¬ 
low 300 mK), which are similar to those used by Her¬ 
schel SPIRE (Griffin et al. 2002, 2003). 

The use of dichroic filters allows BLASTPol to ob¬ 
serve simultaneously in three frequency bands centered 
at 250, 350, and 500 pm. Unlike ground based tele¬ 
scopes it is not restricted to observe through narrow 
windows in the atmospheric transmission spectrum. In¬ 
stead, BLASTPol observes in three wide frequency bands 
(A/// ~ 30%), which bracket the peak of 10—20K 
thermal dust emission. A metal mesh polarizing grid is 
mounted in front of each detector array. The polarizing 
grid is patterned so that each adjacent detector samples 
only the vertical or horizontal component of the incoming 
radiation. In this way a single Stokes parameter {Q or U) 
can be measured in the time it takes light from a source to 
move from one bolometer to an adjacent bolometer (<1 
second for typical scan speeds). A sapphire achromatic 
half-wave plate (hereafter HWP) provides additional po¬ 
larization modulation (Moncelsi et al. 2014). A detailed 
description of the instrument and summary of the ob¬ 
servations will appear in a forthcoming publication (F. 
Angile et al. 2016, in preparation).^^ The present paper 
refers only to observations of Vela G and surrounding 
regions made during the 2012 flight. 

BLASTPol was launched on 26 December 2012. The 
payload rose to an average altitude of 38.5 km above 
sea level and began science operations, taking data un¬ 
til cryogens were depleted 12 days and 12 hours after 
launch. Our selection of target molecular clouds was 
informed by target distance, visibility from Antarctica, 
and cloud brightness. The nearby GMG VelaG was our 
highest priority target. The observations discussed in 
this paper include two types of scans as shown in Figure 
1 (cyan lines). Most of the integration time (43 hours) 
was used to map a “deep” (3.1 deg^) quadrilateral region 
covering four of the hve cloud subregions defined by Hill 
et al. (2011).^“* A further 11 hours were spent mapping a 

See also Matthews et al. (2014) for a description of the 2010 
BLASTPol flight. 

The North region, as defined by Hill et al. (2011), has a sig¬ 
nificant spatial offset from the other four subregions and so was 


larger (~10deg^) area that includes significant regions of 
low dust column where Ay ~ 1 according to extinction 
maps from Dobashi et al. (2005). The larger region was 
observed to reconstruct the map zero-intensity levels. 

Observations were made in sets of four raster scans, 
where the HWP was rotated to one of four angles (0°, 
22?5, 45°, 67?5) after every completed raster scan. A 
complete set of four scans typically required one hour. 
Scans were made while the source was rising and setting 
to maximize the range of parallactic angle. 

As discussed in Section 3.2, the telescope beam was 
much larger than predicted by our optics model, with 
signihcant non-Gaussian structure. We smoothed our 
data to achieve an approximately round beam having a 
FWHM of 2(5 at 500 ^m. 

3. DATA ANALYSIS 

In this section we give a brief overview of the BLAST¬ 
Pol data analysis pipeline and highlight differences from 
the 2010 data pipeline described in Matthews et al. 
(2014). An in-depth description of the data reduction 
pipeline and iterative mapmaker TOAST will be given in 
a forthcoming paper (S. Benton et al. 2016, in prepara¬ 
tion). 

3.1. Time Domain Preprocessing 

Standard techniques were applied to the bolometer 
time-ordered data (hereafter TOD) to remove detec¬ 
tor spikes (mostly due to cosmic rays), deconvolve 
the bolometer time constant, and remove an elevation- 
dependent feature (see Matthews et al. 2014). The data 
were further preprocessed by fitting and removing an ex¬ 
ponential function fit to each detector’s TOD in the first 
30 seconds after a HWP rotation or a telescope slew. A 
high-pass filter with power-law cutoff was used to whiten 
noise in the TOD below 5 mHz. Temporal gain variations 
were removed using the DG voltage level of each detec¬ 
tor and periodic measurements of an internal calibration 
source. Pixel-to-pixel detector gain variations were cor¬ 
rected by frequent observations of the bright compact 
source IRAS 08470-4243. 

Telescope attitude was reconstructed using pointing 
solutions generated from the BLASTPol optical star 
camera,^® with payload rotational velocities from gy¬ 
roscopes used to interpolate between pointing solutions 
(Pascale et al. 2008). Data having pointing uncertainties 
> 5" were discarded. The final on-sky pointing solution 
was calibrated to match the astrometry of publicly avail¬ 
able Herschel SPIRE maps^® at the same wavelength. 

3.2. Beam Analysis 

The BLASTPol 2012 beam differs from the beam pre¬ 
dicted by our optics model. It has multiple elongated 
peaks, and the relative power in each peak varies from 
detector to detector. BLASTPol filters were designed for 
near-space conditions and the telescope far field is sev¬ 
eral kilometers away so it was not possible to map the 

not included in our deep scan region. 

BLASTPol flew two redundant star-boresight optical star 
cameras during the 2012 flight, but one experienced a harddrive 
failure six hours after the launch (Galitzki et al. 2014). 

http://www.cosmos.esa.int/web/herschel/ 

science-archive 
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Figure 1. BLASTPol 500/im I map of VelaC and surroundings. Cyan lines show the boundaries of the two raster scan regions used to 
make the maps in this paper: the region marked with cyan dashes was observed for 43 hours, while the solid cyan lines show a larger region 
covering VelaC and surrounding diffuse emission, which was observed for 11 hours in total. Also shown are the locations of the regions 
used in the diffuse emission subtraction for the BLASTPol /, Q, and U maps as described in Section 3.5. The region labeled C is used 
in the “conservative” diffuse emission subtraction method. The “aggressive” method used the two regions labelled Al and A2 . Cloud 
subregions defined by Hill et al. (2011) are indicated in white contours. The region outlined in blue shows the “validity region” where the 
null tests discussed in Section 3.6 were passed, and where both diffuse emission subtraction methods discussed in Section 3.5 are valid; 
only polarization measurements within this validity region are used for science analysis in later sections. The red circle shows the area near 
RCW 36 excluded from our polarization analysis because of large Stokes 7, Q, and U null test residuals. 


far-field beam at sea level (Galitzki et al. 2014). Instead 
the beam had to be inferred from in-flight measurements 
of astronomical objects. 

Our 2012 instrument beam model was defined in tele¬ 
scope coordinates and was informed by observations of 
two objects: IRAS 08470—4243, a warm compact dust 
source located in the Vela Molecular Ridge; and limited 
observations of the planet Saturn made on 27 December 
2015. IRAS 08470—4243 was observed every 4—8 hours, 
with reasonable coverage for all detectors, but it is not 
a point source at BLASTPol resolution. BLAST obser¬ 
vations of IRAS 08470-4243 in 2006 found a FWHM of 
^40" (Netterfield et al. 2009). Saturn has a radius of 
6 X 10"* km, which corresponds to an angular size of <20", 
considerably smaller than the BLASTPol 2012 beam. 
Saturn was only observed early in the flight at telescope 
elevations of <25° and was only fully mapped by the 
bolometers near the center of the focal plane arrays. 
Three elliptical Gaussians were fit to the Stokes I maps 
of IRAS 08470—4243 and Saturn. The free parameters 
were the Gaussian amplitudes, centroids, widths, and po¬ 
sition angles. Only pixels above an intensity threshold 
of 20% of the peak intensity for IRAS 08470—4243 and 
7.5% of the peak intensity for Saturn were used in the 
fits. The final 2012 instrument beam model used the 
centroid positions, amplitudes, and position angles from 
the fits to IRAS 08470—4243 and the Gaussian widths 
calculated from the fits to Saturn. 


Next, an on-sky beam model was computed for the 
VelaG observations. The on-sky beam model is a time- 
weighted average of the instrument beam model rotated 
by the angle between the telescope vertical direction and 
Galactic north for each raster scan. The resulting aver¬ 
age beam model for Vela G is shown in the left panel of 
Figure 2. A single elliptical Gaussian fit to this beam 
model gives FWHMs of 130" by 64", at position angle 
—51°. This beam is significantly larger than the expected 
diffraction limit of the telescope (FWHM=60"). 

Lucy-Richardson (LR) deconvolution was previously 
used to correct a larger-than-expected beam from the 
BLAST 2005 flight (Roy et al. 2011). Here we used an 
iterative LR method and our beam model to deconvolve 
a simulated map consisting of a single Gaussian source 
with FWHM = 145". The deconvolved map from this 
step when applied as a smoothing kernel to convolve the 
original beam model should restore the 145" Gaussian. 
The success of this step can be judged from the right 
panel in Figure 2: it does produce a single-peaked source 
that is approximately round (FWHM=144" by 157"). 
This same smoothing kernel was used to convolve the /, 
Q, and U maps of Vela G, with a resulting resolution of 
approximately 2^5. 

3.3. Instrumental Polarization 

To determine the instrumental polarization (the po¬ 
larization signal introduced by the instrument hereafter 
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Figure 2. Left: BLASTPol 500/i.m beam model for the VelaC map. Right: BLASTPol 500 ^m beam model for the VelaC map after 
convolution with the smoothing kernel discussed in Section 3.2. Contour levels (cyan) are 25, 50, and 75 % of the peak brightness. The 
dashed blue lines in each image show the FWHM of a fit to an elliptical Gaussian. FWHMs of the fitted Gaussians are 130'^ by 64" for 
the VelaC beam model (left panel) and 144" by 157" for the smoothed beam model (right panel). 



referred to as IP) we followed methods described in 
Matthews et al. (2014). In brief, the set of observations 
of Vela C was split into two bins based on parallactic 
angle, and maps were produced for each detector indi¬ 
vidually using the naivepol mapmaker (Moncelsi et al. 
2014). The measured polarization is a superposition of 
one component fixed with respect to the sky and the IP, 
which is fixed with respect to the telescope. By compar¬ 
ing the polarization measurements at different parallactic 
angles the IP of each bolometer could be reconstructed. 
These IP terms were then removed during the mapmak¬ 
ing stage (Section 3.4). The 500/rm array has an aver¬ 
age IP amplitude of 0.53%. To check the effectiveness of 
the IP correction the VelaC data were divided into two 
halves and IP estimates derived from the first half of the 
data were used to correct the second half of the data. 
By measuring the IP of the “corrected” second half of 
the VelaC data we estimate that the minimum value of 
the fractional polarization p measurable by BLASTPol at 
500 ^m is 0.1%. 

3.4. TOAST 

Maps were made using TOAST (Time Ordered As¬ 
trophysics Scalable Tools),a collection of serial and 
OpenMP/MPI parallel tools for simulation and map 
making. Specifically, the generalized least-squares solver 
was used, which iteratively inverts the map-maker equa¬ 
tion using the preconditioned conjugate gradient method 
(see Cantalupo et al. 2010, a predecessor to TOAST). 
The map-maker’s noise model was estimated using power 
spectra from observations of faint dust emission in the 
constellation Puppis (map centered at I = 239.0°, 
b = —1.7°), with simulated astrophysical signal sub¬ 
tracted. The noise model is consistent with white noise 
plus (1//)“ correlations that level out at low frequency 
due to data preprocessing. Correlations between detec- 

http://tskisner.github.io/TOAST 


tors and non-stationarity of the noise were not required 
by the model. Instrumental polarization (Section 3.3) 
was removed as per Matthews et al. (2014). Per-pixel 
covariance matrices for Stokes /, Q, and U were esti¬ 
mated as the 3x3 diagonal block of the full pixel-pixel 
covariance matrix. Noise-only maps, both simulations 
and null tests (see Section 3.6), are consistent with the 
estimated covariances, up to a constant scaling factor due 
to unmodeled noise. A pixel size of 10" was used for all 
signal and covariance maps. 

3.5. Dijfuse Emission Subtraction 

To study the polarization properties and magnetic field 
morphology of Vela C it is necessary to isolate polarized 
dust emission originating in VelaC from the diffuse po¬ 
larized emission associated with Galactic foreground and 
background dust as well as the Vela Molecular Ridge. 
This separation requires extra care as previous studies 
show that diffuse sightlines, which may be used to es¬ 
timate the foreground/background polarized emission, 
tend to have a higher average polarization fraction than 
dense cloud sightlines. In particular, Planck observations 
show that in the more diffuse clouds there is a range of 
p values with the maximum reaching 15 to 20%, while 
such high values are not seen in the higher column den¬ 
sity clouds (e.g., Orion, Ophiuchus, Taurus), where the 
average p is consistently lower (Planck Collaboration Int. 
XX 2015). Polarized emission from diffuse dust along 
dense cloud sightlines could therefore contribute signifi¬ 
cantly to the overall polarization measured. In this sec¬ 
tion we present two different diffuse emission subtraction 
methods, one conservative and one more aggressive with 
respect to diffuse emission removal. 

In the conservative method for diffuse emission sub¬ 
traction, we considered most of the emission surrounding 

It should also be noted that observations made by BLASTPol 
are inherently differential measurements, and thus the map zero- 
intensity level is uncertain. 
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Vela C as defined by Hill et al. (2011) to be associated 
with the cloud. The Hill et al. (2011) Vela C cloud sub- 
regions are overplotted (solid white lines) on a map of 
500 total intensity in Figure 1. The zero-point for 
the intensity of the cloud emission was set by specify¬ 
ing a region with relatively low intensity near Vela C 
and assuming that emission in this reference region also 
contributes to sightlines on the cloud with spatial unifor¬ 
mity. This low flux region is labeled “C” in Figure 1. We 
calculated the average Stokes I, Q, and U in that region, 
and the appropriate mean flux was then subtracted from 
each of the maps. The result was a set of maps of Vela C 
emission isolated from the Galaxy, assuming a minimal, 
uniform contribution from foreground and background 
emission. 

In the aggressive method we considered most of the dif¬ 
fuse emission surrounding Vela C to be unassociated with 
the cloud, and accordingly a higher level of flux needed to 
be removed to isolate the cloud. Furthermore, we noted 
that there is significantly more emission to the south of 
Vela C than to the north; thus it is reasonable to assume 
that the true map of the region consists of the Vela C 
cloud superimposed on a varying Galactic emission pro¬ 
file. In this method of referencing, we defined two regions 
closely surrounding the cloud (marked “Al” and “A2” in 
Figure 1) and performed 2-dimensional linear-plane fits 
to the Stokes /, Q, and U maps excluding all map pixels 
except those located in regions Al and A2. The three 
free parameters in these fits were the linear slopes of the 
plane in the directions tangent to I and b and a map off¬ 
set. The equations for each of the resulting plane-fits to 
the I, Q, and U maps were used to specify the intensity to 
be subtracted from each pixel in the maps. Note that for 
regions far from Vela C, the linear approximation of the 
Galactic emission profile breaks down, leading to an in¬ 
appropriate extrapolation. Therefore we defined an area 
within which the linear fit referencing method is valid 
(blue quadrilateral in Figure 1), bounded on the north 
and south by the reference regions Al and A2, and on 
the east and west by the edges of the well-sampled por¬ 
tion of the map . This “validity” area roughly coincides 
with the four southernmost regions of Hill et al. (2011). 
We note that some of the emission in Al and A2 might 
in fact be associated with VelaC, so this method is likely 
to over-subtract the diffuse dust emission. 

The true I, Q, and U maps of Vela C probably ex¬ 
ist somewhere between our most extreme physically rea¬ 
sonable assumptions corresponding to the conservative 
and aggressive diffuse emission subtraction methods. In 
this paper we present results for an “intermediate” dif¬ 
fuse emission subtraction method, derived by taking the 
arithmetic mean of the /, Q, and U maps from the ag¬ 
gressive and conservative methods. Most of our analyses 
are then repeated using the aggressive and conservative 
methods as a gauge of the uncertainties associated with 
the diffuse emission subtraction. 


3.6. Null Tests 

To characterize possible systematic errors in our data, 
we performed a series of null tests, which are described in 
detail in Appendix A. In these, we split the 250 fim obser¬ 


vations^® into two mutually exclusive sets. If the polar¬ 
ization parameters from the two independent data sets 
agree, we can conclude that the impact of systematics 
is small, and the uncertainties are properly character¬ 
ized by Gaussian errors produced by TOAST. The four 
methods of splitting are: data from the left half of the 
array vs. the right half; data from the top half of the 
array vs. the bottom half; data from earlier in the flight 
vs. later in the flight; and alternating every other scan 
set sequentially throughout the flight. 

For each null test we made separate maps of I, Q, and 
U, which were then used to calculate residual maps of the 
polarized intensity P, the polarization fraction p, and the 
polarization-angle ip as described in Appendix A. (The 
quantities P, p, and ip are defined in Appendix B). If 
our data had no systematic errors we would expect to 
see uncorrelated noise in the residual maps. For P and 
p if the residuals were less than one-third of the signal 
in a given map pixel then that pixel was said to pass the 
null test. For ip the residuals had to be less than 10° to 
pass the null test. We examined each of the four null 
tests listed above for each of the two methods of diffuse 
emission subtraction (Section 3.5) in P, p and ip giving a 
total of 24 checks that our measured polarization signal 
is significantly above the systematic uncertainty level. 

We found that our map passed these tests for the ma¬ 
jority of sightlines inside the cloud region shown in Figure 
1 (blue solid line). The exceptions occurred in regions 
where the fractional polarization was small, so that a 
comparison of the scale of the polarization signal to the 
scale set by residuals in the null tests resulted in an ap¬ 
parent failure. The fact that we saw null test failures 
correlating with low p, but not with absolute difference 
in the null test / maps led us to the interpretation that 
the apparent low signal level compared to the null test 
residual is due to decreased signal and not increased sys¬ 
tematic uncertainties. We did see significant structure 
in the null test residual maps of Q and U near the com¬ 
pact H II region RGW 36, which coincided with null test 
residuals of one-fourth p, though the residuals in ip were 
much smaller than 10°. These p measurements techni¬ 
cally pass the null test criteria, but the systematic errors 
are larger than the statistical errors. We conclude that 
for the validity region shown in Figure 1 the null tests 
are passed, with the exception of a circular area centered 
on RGW36 (Z =265.15°, b = 1.42°within a radius of 4'). 

4. BLAST-POL POLARIZATION MAPS 

In this section we present maps of the Stokes param¬ 
eters I, Q, and U, linearly polarized intensity (P), and 
polarization fraction (p = P/I). The polarization de¬ 
scriptors and covariances used in our analysis are sum¬ 
marized in Appendix B. We also present maps of $, the 
inferred orientation of magnetic field projected onto the 
plane of the sky, which is assumed to be the orienta¬ 
tion of the polarization of the dust emission (described 
by Ip) rotated by 90°, and the localized dispersion in the 
polarization-angle (S'). 

Because p and P are constrained to be positive any 
noise in the Q and U maps will tend to increase the 

As discussed in Appendix A the BLASTPol 250 pm observa¬ 
tions of Vela C are better suited than the 500 pm data for perform¬ 
ing null tests. 
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Figure 3. BLASTPol 500/im VelaC maps of 7, Q, U, and the total polarized intensity P (which is debiased as described in Section 
4). The color scale units are MJysr“^. Contours indicate I levels of 46, 94, 142, and 190 MJysr“^, and the gray outlines indicate cloud 
subregions covered by our observations as in Figure 1. 


measured polarization. Accordingly, we crudely debias 
p and P according to 


Pdb = 

\/p" 

- 

(1) 

Pdb — 

JP^ 

- o-|,, 

(2) 


(Wardle & Kronberg 1974). This method of debiasing 
is appropriate only where Up is small compared with p 
(Montier et al. 2015). We note that the median value of 
p/<Jp in our map is ~ 25, so for most of our map this 
debiasing method is applicable. 

4.1. Diffuse Emission Subtracted Maps of I, Q, and U, 
and Derived Maps of P, and p 

Figure 3 shows Vela C 500 pm maps for the three Stokes 
parameters /, Q, and U. The maps have been smoothed 
to 215 resolution, as described in Section 3.2, and use the 
intermediate diffuse emission subtraction method (Sec¬ 
tion 3.5). Overlaid in gray are the outlines of the subre¬ 


gions of Vela C as defined in Hill et al. (2011) and labelled 
in Figure 8. The BLASTPol I map peaks at the location 
of RCW36. 

Also included in Figure 3 is the derived map of the po¬ 
larized intensity (P), which generally shows some signal 
where there is cloud emission. However, the correspon¬ 
dence is certainly not perfect, and varies considerably 
across the map. For example, along most of the Centre- 
Ridge there is a corresponding peak in the P map along 
the main ridge. In the South-Ridge there are peaks at 
similar locations in the P and I maps. But in the South- 
Nest, prominent areas of polarized emission are only seen 
around the edge of the cloud structure seen in /. There 
are also some regions of significant P that stand out less 
in /, for example, along the north edge of the Centre- 
Ridge. 

Figure 4 shows the polarization fraction {p = P/I) 
for each of the three different diffuse emission subtrac¬ 
tion choices discussed in Section 3.5. The conservative 
diffuse emission subtraction (top panel) results in p that 
is lower on average than from the aggressive diffuse emis- 
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Figure 4. BLASTPol 500/im maps of p obtained using different 
methods for separating the polarized emission of VelaC from that 
of diffuse background/foreground dust (Section 3.5): conservative 
method (top panel); aggressive method (middle panel), and inter¬ 
mediate method (bottom panel). Only sightlines where p > 3crp 
and 7 > 0 are shown. The p maps shown have been debiased 
using the methods described in Section 4. Gray contours indicate 
I levels of 46, 94, 142, and 190 MJysr“^, and the white outlines 
indicate the four VelaC cloud subregions as in Figures 1 and 8. 


sion subtraction (middle panel). This is expected as p is 
commonly observed to increase for regions of low dust 
emission. Thus, compared to the aggressive subtraction 
method that uses regions closer to the cloud with lower 
average p, the conservative method removes more P rel¬ 
ative to I. The bottom panel shows the p map result¬ 
ing from the intermediate diffuse emission subtraction 
method. Unless otherwise specified, for the remainder of 
the paper p, ip, and are calculated from Stokes param¬ 
eter maps using the intermediate diffuse emission sub¬ 
traction method. 

The mean value of p in our map is 6.0% with a me¬ 
dian of 3.4 %. For the map pixels within the dense cloud 
subregions defined by Hill et al. (2011) the mean po¬ 
larization fraction is 3.5% with a median of 3.0% and 
a standard deviation of 2.4%. Previous submillimeter 
polarization maps having spatial coverage correspond¬ 
ing to the scales of entire clouds have yielded roughly 
similar values. Specifically, after subtracting the back¬ 
ground/foreground emission, Planck Collaboration Int. 
XXXIII (2016) found mean 850 pm polarization frac¬ 
tions of 1.8%, 5.0%, and 6.1% for three nearby molecular 
clouds, while 450 pm polarization maps of four GMCs 
made by Li et al. (2006) with SPARO at the South Pole 
yield a mean polarization fraction of 2.0%. Our p map 
shows behavior that is broadly consistent with expecta¬ 
tions from the P map. Values of p tend to decrease with 
increasing /, but there is not a one-to-one anticorrelation 
between p and /. 

4.2. Inferred Magnetic Field Direction 

Figure 5 shows a detailed view of the magnetic field 
orientation projected onto the plane of the sky $, as 
inferred from the BLASTPol 500 pm data. This figure 
uses a “drapery” pattern produced using the line inte¬ 
gral convolution method of Cabral & Leedom (1993) su¬ 
perimposed on the BLASTPol 500 pm I map.^*^ Dotson 
(1996) showed that there is significant ambiguity in infer¬ 
ring the magnetic field lines from polarization data, par¬ 
ticularly as polarization maps can sample multiple cloud 
structures along the line of sight, each potentially having 
a different magnetic field orientation. The drapery im¬ 
age is presented solely to show the range of orientations 
of 4), and to give a sense of the range of spatial scales 
probed by BLASTPol. Figure 6 shows 4> as a series of 
line segments (approximately one line segment per 215 
BLASTPol beam). 

The projected cloud magnetic field direction appears 
to change across VelaC: at low Galactic latitudes the 
field is mostly perpendicular to the main cloud elongation 
direction, while at higher Galactic latitudes it bends to 
run mostly parallel to the cloud elongation direction. We 
also see some sharp changes in <I>, most noticeably in the 
South Nest, and near the compact H II region RGW 36. 

Figure 7 shows that the dispersion in magnetic field 
orientation across the Vela G cloud is 28°. Novak et al. 
(2009) calculated dispersions on similar spatial scales by 
combining the large-scale GMG polarization maps of Li 
et al. (2006) with higher angular resolution submillime¬ 
ter polarimetry data. They obtained 27°-28°, nearly 

This visualization is produced with the same code used in 
Planck Collaboration Int. XXXV (2016) and will be further dis¬ 
cussed in a forthcoming paper by D. Falceta-Gongalves et al. 
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Figure 5. BLASTPol 500 /im I map with the inferred plane of the sky magnetic field component (<I>) overlaid as a “drapery” image (only 
regions where I > 0 are shown). The drapery pattern is produced using the line integral convolution method (Cabral & Leedom 1993) 
and indicates the orientation of the magnetic field as projected on the plane of the sky. Note that this drapery pattern was made from all 
of the $ data with no masking of sightlines having large uncertainties in $. This image is meant show the level of detail available in the 
BLASTPol $ maps, but should not be used for quantitative analysis. 


the same result. In future publications we will present 
statistical studies of the correlations between magnetic 
field orientation, filamentary structure, and cloud veloc¬ 
ity structure. 

4.3. Polarization Angle Dispersion Function 

To quantify the disorder of $ in our VelaC maps at 
small scales we calculate the polarization-angle disper¬ 
sion function S', implementing the formalism described 
in Section 3.3 of Planck Collaboration Int. XIX (2015). 
For each pixel in our map, S is defined as the rms devia¬ 
tion of the polarization-angle '0 (x) for a series of points 
on an annulus of radius 6: 

1 ^ 

= ( 3 ) 

i=0 

where S is the length scale of the dispersion, x is the po¬ 
sition for which we evaluate the polarization-angle dis¬ 
persion and 

Sxi =-Ip (x) -'ll) (^x + . (4) 

Because S is always positive it is biased due to noise. 
We debias using the standard formula 

Sj,iS) = S^S) - al (5) 

where cr| is the variance of S. 


Figure 6 shows S for S = 2'.5 (^0.5pc), the smallest 
scale that can be resolved with our smoothed beam. 
(Hereafter we refer to S (2(5) as S). The most striking 
features in the S map correspond to regions of sharp 
changes in $, which is indicated with line segments. 
These high dispersion regions sometimes occur near the 
locations of dense filaments (for example, the sharp bend 
in the South Nest). More often they correspond to sight¬ 
lines of lower than average p and do not appear to be 
coincident with any prominent cloud feature in I. 

4.4. Polarization Map Sampling and Sightline Selection 

Our BLASTPol polarization maps were calculated 
from Stokes parameters smoothed to a resolution of ^2(5 
as discussed in Section 3.2. The resulting polarization 
maps were then sampled every 70" to ensure at least 
Nyquist sampling. In total there are 4708 projected mag¬ 
netic field sightlines over the validity region defined in 
Section 3.5. 

In the following sections we attempt to model the po¬ 
larization fraction p as a function of N, T, and S. For 
these detailed studies we restrict our analysis to sight¬ 
lines that encompass the dense cloud regions as defined 
by Hill et al. (2011). These sightlines are better probes of 
the polarization structure in the cloud and are less sen¬ 
sitive to systematic uncertainties in our ability to sepa¬ 
rate the polarized emission emitted by diffuse dust fore¬ 
grounds/backgrounds from the polarized emission emit¬ 
ted by dust grains in VelaC. 
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Figure 6. BLASTPol 500 fim map of the dispersion in the polarization-angle (S) in degrees on 0.5 pc scales (<5 =2.5') as defined in Section 
4.3. The S map is shown where S > 3crs. Line segments show the orientation of the magnetic field as projected on the plane of the sky 
derived from the BLASTPol 500/rm data. The <I> measurements are shown approximately every 2f5. Contours indicate 500/^m I intensity 
levels of 46, 94, 142, and 190 MJysr“^. 


To ensure a robust sample, we use only p values that 
are large enough to be unaffected by uncertainties in in¬ 
strumental polarization removal {p > 0.1%, see Section 
3.3), and for which we have at least a Scr detection of 
polarization (p > 3ap), which corresponds to an un¬ 
certainty in the polarization angle < 10°. To en¬ 
sure that the polarization values are not dependent on 
our choice of diffuse emission subtraction method we re¬ 
quire that pint > 3|pint -PconI and Pint > 3|pi„t -Paggi, 
where Pconj Pagg, and pint are the polarization fraction 
values calculated using the conservative, aggressive, and 
intermediate diffuse emission subtraction methods re¬ 
spectively (see Section 3.5). Similarly, we require that 
li^int-'0con| < 10° and |i/’int-i/’aggi < 10°. We also ex¬ 
clude sightlines from a 4' radius region near RCW 36 as 
these show residual structure in our null tests (see Sec¬ 
tion 3.6). In total 2488 out of a 3056 possible Hill et al. 
(2011) sightlines meet these criteria. For our analysis of 
p vs. N, T, and S in Sections 6 and 7 we also require at 
least 3-cr measurements of N, T, and S where the errors 
on N and T are derived from the SED fit covariance ma¬ 
trices (see Section 5). This results in a final sample of 


2378 sightlines. 

5. COLUMN DENSITY AND TEMPERATURE MAPS 
DERIVED FROM Herschel SPIRE DATA 

To derive column density and dust temperature 
maps we used publicly available Herschel SPIRE and 
PACS data. SPIRE uses nearly identical filters to 
BLASTPol, but has higher spatial resolution (FWHM of 
17'.'6, 23"9, and 35"2 for the 250, 350, and 500 pm bands, 
respectively). Data taken with the PACS instrument in 
a band centered at 160 pm (FWHM of 13'.'6) were used 
to provide additional sensitivity to warm dust. Her¬ 
schel maps were generated using Scanamorphos (Rous¬ 
sel 2013) and additional reduction and manipulation was 
performed in the Herschel Interactive Processing Envi¬ 
ronment (HIPE version 11) including the Zero Point Cor¬ 
rection function for the SPIRE maps. The resulting Her¬ 
schel maps were smoothed to 35'.'2 resolution by convolv¬ 
ing with Gaussian kernels of an appropriate size and then 
regridding to match the 500 pm map. 

Similar to the diffuse emission subtraction described in 
Section 3.5, we attempted to separate the Galactic fore¬ 
ground and background dust emission from the emission 
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Figure 7. Histograms of the BLASTPol 500 inferred magnetic 
field direction $ for all VelaC sightlines (red) and sightlines lying 
inside the Hill et al. (2011) subregions (blue). Sightlines included in 
these histograms have cr$ < 10° and both ji/jint —V’conI < 10° and 
I'^Aint — '0agg| < 10° (see Section 4.4). The standard deviation of 
each distribution is given at upper-left. 

of VelaC. As the regions used to define the diffuse emis¬ 
sion subtraction in Section 3.5 were not covered by the 
Herschel map we defined four alternate “diffuse emission 
regions” (see Figure 8 top panel). These regions were 
presumed to contain little emission from dust in Vela C 
and thus they are reasonable representations of the con¬ 
tribution due to diffuse dust emission. For the initial 
analysis described below, the mean intensity in Region 
1 was subtracted from each of the 160, 250, 350, and 
500 ^m maps, and the maps were then further smoothed 
to match the 2l5 resolution of the BLASTPol maps. 

Modified blackbody SED fits were made for each map 
pixel using the methods described in Hill et al. (2009, 
2010, 2011) and using the dust opacity law of Hildebrand 
(1983) with a dust spectral index (3 = 2. The resulting 
column density (N) and dust temperature (T) maps are 
shown in Figure 8 (middle and bottom panels, respec¬ 
tively). It should be noted that above a temperature of 
^20 K, the dust emission is expected to peak at wave¬ 
lengths shorter than 160 /im. For these warmest sight¬ 
lines our estimates will have a higher degree of uncer¬ 
tainty. The derived N and T maps were visually com¬ 
pared to the higher resolution column density and tem¬ 
perature maps from Hill et al. (2011), which did not in¬ 
clude a diffuse emission subtraction. Our maps are in 
close agreement with the Hill et al. (2011) maps for col¬ 
umn density sightlines where VelaC emission is strong 
compared to the diffuse emission component. Note that 
we computed maps of the column density of hydrogen nu¬ 
clei while Hill et al. (2011) calculated the column density 
of H 2 . 

Much of the analysis in the present paper focuses on 
comparisons between parameters such as polarization 
fraction p, N, and T. From Figure 8 we see that N and 
T are strongly anti-correlated. Similar trends were noted 
by Palmeirim et al. (2013) in their Herschel study of a 
cold cloud in Taurus. We interpret this as a result of radi¬ 
ation shielding in the densest parts of the cloud. This in¬ 
terpretation can be tested by examining a plot of 250 pm 
intensity vs. 500 pm intensity, as shown in Figure 9. In 


Figure 8. Herschel VelaC 500 fim intensity (top panel, 
FWHM = 35'.'2), column density [N, middle panel, FWHM = 2l5), 
and dust temperature (T, bottom panel, FWHM = 2l5). The 
N and T maps were derived from Herschel data using the meth¬ 
ods described in Section 5. Numbered quadrilaterals correspond 
to different diffuse emission regions for which the average intensity 
is indicated in Figure 9. Note that the mean intensity in Region 1 
was subtracted from each of the 160, 250, 350, and 500 ^m maps 
before SED fitting. The solid black polygons (labeled in the top 
panel) correspond to the cloud subregions as defined in Hill et al. 
( 2011 ). From left to right these are: the South Nest, a region 
of many overlapping filaments; the South Ridge, dominated by a 
single dense filament; the Centre Nest; and Centre Ridge, which 
contains the ionizing source(s) powering the compact H II region 
associated with RCW 36. Hill et al. (2011) also include an addi¬ 
tional region, designated North, that was not covered in the deep 
BLASTPol survey of VelaC. 

this figure there is a noticeable bend in the otherwise lin¬ 
ear relationship between the two intensities. Since sub¬ 
millimeter dust emission in molecular clouds is typically 
optically-thin, larger intensity at either wavelength cor¬ 
responds to higher column density. However, beyond the 
bend we notice that the slope of the 500 pm intensity vs. 
250 pm intensity relation decreases. The simplest expla¬ 
nation is that the dust in denser regions of the cloud is 
colder, due to radiation shielding. 

An alternative interpretation of the bend seen in Fig¬ 
ure 9 is to hypothesize a uniformly cold cloud spatially 
superimposed on diffuse emission from warmer dust. To 
explore this possibility we examined the location of each 
diffuse region on Figure 9 relative to the bend in the ob¬ 
served curve of 250 pm vs. 500 pm intensity. Subtracting 
the diffuse emission flux essentially sets a new origin for 
this graph and is equivalent to the diffuse emission sub¬ 
traction discussed in Section 3.5, leaving only emission 
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Figure 9. Median values of Herschel 250//m intensity in bins of 
Herschel 500 intensity for the South Nest region in Vela C (as 
labeled in Figure 8). The error bars correspond to the standard 
deviation of the 250 pm intensity values in each bin. Black lines 
correspond to the expected intensity ratios for uniform temperature 
dust. Diamonds indicate the average 250 pm and 500 pm intensities 
for the four numbered diffuse emission regions indicated in the top 
panel of Figure 8 (from left to right these indicate regions 1, 2, 3, 
and 4). Error bars show the standard deviation of intensity values 
in each region. 

from dust grains in the Vela C cloud. As can be seen 
from Figure 9, the diffuse emission regions, even very ag¬ 
gressively placed ones, reposition the origin to locations 
significantly below the bend in the curve, indicating that 
the observed T and N anticorrelation is intrinsic to the 
Vela C molecular cloud. As a further check the SED fits 
described above were redone using diffuse emission Re¬ 
gions 2,3, and 4 as the reference regions, instead of using 
Region 1 (see Figure 8). The corresponding N maps are 
very similar to the one shown in Figure 8, especially for 
the densest regions. 

6. DEPENDENCE OF POLARIZATION FRACTION ON 
N AND T 

Before considering the polarization fraction p, we first 
attempt to separate sightlines that show significant heat¬ 
ing from sources internal to VelaC from sightlines that 
appear to be predominantly heated by the interstellar ra¬ 
diation field (ISRF). The polarization properties of sight¬ 
lines near a source of intense radiation, such as the com¬ 
pact H II region RCW36 in VelaC, might differ from the 
polarization properties of cloud sightlines where star for¬ 
mation is at an earlier stage. The presence of a bright 
radiation source might affect the efficiency of radiative 
torques in aligning dust grains with respect to the local 
magnetic field. Also, the presence of expanding ionized 
gas in H II regions can alter the magnetic field geometry, 
for example as seen in SPARO observations of the Carina 
Nebula (Li et al. 2006). 

Figure 10 shows T vs. logV for sightlines selected as 
discussed in Section 4.4. (Note that throughout this pa¬ 
per log refers to log^o-) As discussed in Section 5 the 
ISRF can more easily penetrate sightlines of low column 
and therefore average temperatures of low N sightlines 
tend to be higher. Figure 10 generally shows decreasing 



Figure 10. Dust temperature (T) vs. logA^ for all Vela C sight¬ 
lines in the subregions defined by Hill et al. (2011). Blue diamonds 
show sightlines that were rejected by an iterative application of 
Chauvenet’s criterion. These 143 sightlines appear to be heated 
by the compact H II region RCW 36. The other 2235 sightlines 
(crosses) appear to be heated only by the interstellar radiation field 
(ISRF). The red line corresponds a fit to all ISRF-heated sightlines, 
as described in Section 6. 

T with increasing log N, however it also shows that a mi¬ 
nority of sightlines have temperatures lying well above 
this approximately linear trend. We fit the equation 
T = alogN + b, using Chauvenet’s criterion (Chau- 
venet 1863) iteratively to remove outliers (diamonds in 
Figure 10). The 143 sightlines rejected as outliers are lo¬ 
cated near the compact H II region RCW 36 (upper panel 
of Figure 11). These sightlines appear to be heated by 
the H II region, yielding temperatures lying above the 
trend seen for ISRF heated sightlines. 

Figure 12 shows the dependence of p on V and T for 
ISRF-heated sightlines (left side, top and middle pan¬ 
els respectively). In general p decreases with increasing 
N and increases with increasing T. To quantify the de¬ 
pendence of p on V we fit a model of the form 

p = CV"" (6) 

where, C and are constants. This is equivalent to a 
linear fit in logarithmic space 

logp = UN logV -b C. (7) 

Via a fit to Equation (7), we find that aAr=—0.58 ± 0.02. 
Each measurement of logp is given equal weight in our fit. 
By giving each data point equal weight (equal fractional 
error in p) we are assuming that the deviations of the 
logp data points from the fit described in Equation (7) 
are caused by additional dependences of p on other quan¬ 
tities, rather than uncertainties in our measurements of 
logp. This assumption is reasonable, because our po¬ 
larization measurement uncertainties are generally quite 
small. For example, the median signal-to-noise of our 
p measurements is 36, and the median signal-to-noise of 
the N measurements for these same sightlines is even 
higher. The uncertainties on our fitted parameters are 
calculated using the bootstrapping method with replace¬ 
ment (Press et al. 1992), repeating the fits for each of 
10,000 random selections. The standard deviation of the 
derived power-law exponents is used as an estimate of 
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Figure 11. Color-coded plot of p over the range of 0.002 to 0.100 
for all RCW 36 heated (top panel) and ISRF heated (middle panel) 
BLASTPol sightlines that pass the criteria described in Section 4.4. 
The bottom panel shows which is p for the ISRF heated 

sightlines decorrelated from N and S using Equation (11) in Section 
7.2. If Equation (10) accounted for the entire variation of p, then 
the value of would be constant at 0.029. The background 

image is /soo- 


their uncertainty. 

Similarly, for p vs. T (Figure 12 middle left panel), 
we fit to the relation logp = PtT + ci, and find that 
Pt = 0.125zb 0.005, which implies that p oc exp(0.29r). 
However, Figure 10 shows that N and T are highly cor¬ 
related for ISRF heated sightlines. We can remove the 
correlation of p with N by computing: 


[N] 

Pi = Pi 



( 8 ) 


where , Pi, and Ni are the *th decorrelated p measure¬ 
mentand original p and N measurements, respectively, 
and N is the median value of N for our sightlines. The 
bottom left panel of Figure 12 shows vs. T. By re¬ 
moving the anticorrelation of p with iV, we also remove 
any correlation with T. Thus it appears that there is no 
correlation between p and T that is independent of the 
correlation between p and N. 

For the sightlines that show significant heating 
from RCW 36 we see a similar decrease in p with 
increasing N and find a power-law exponent of 
aN =—0.78±0.06 (Figure 12 top right panel). However, 
for these sightlines there is no apparent correlation be¬ 
tween p and T (Figure 12 middle right panel). 


7. DEPENDENCE OF POLARIZATION FRACTION ON 
N AND S 

In this section our goal is to build an empirical model 
for the dependence of p on iV and the polarization-angle 
dispersion on 2'.5 (0.5pc) scales 5, for an early stage 
star-forming region. Therefore we only consider ISRF- 
heated sightlines. Additionally, we do not include T as 
a parameter of the empirical model as it was shown in 
Section 6 that the p vs. N and p vs. T correlations are 
degenerate. We choose N rather than T as our indepen¬ 
dent variable because the most natural explanation for 
the N vs. T anticorrelation for ISRF-heated sightlines is 
that the density structure of the cloud determines the 
average temperature of the sightlines, rather than T de¬ 
termining N. 


7.1. Individual Correlations among p, N, and S 

Figure 13 shows the median p (color map) for bins of 
S and N for ISRF-heated sightlines. There is a clear 
decrease of p with increasing N and S. Individual corre¬ 
lations are shown in Figures 12 and 14, and the derived 
associated power-law exponents are listed in Table 1. 

Decreasing p with increasing N has been seen in sub¬ 
millimeter polarization observations of many clouds and 
cores (e.g., Matthews et al. 2001; Li et al. 2006). The 
observed decrease in p is often attributed to either can¬ 
celation of polarization signal for high-V sightlines due 
to more disorder in the magnetic field, or to changes in 
grain alignment efficiency within the cloud. These pos¬ 
sible explanations are discussed further in Section 8. 

In Section 4.3, we showed that VelaC has high values 
of S in localized filament-like regions, where there are 
sharp changes in magnetic held direction. High S de¬ 
pends implicitly on spatial changes in the magnetic held 
locally in the map, and any related changes in the mag¬ 
netic held direction within the volume sampled by the 
BLASTPol beam could lead to lower p. The top panel 
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Figure 12. Two-dimensional histograms showing the correlations between p, N, and T for ISRF-heated sightlines (left) and sightlines 
that show evidence of heating from RCW36 (right). The correlations shown are: polarization fraction (p) vs. column density (N) (top 
panels); p vs. dust temperature (T) (middle panels); and the polarization fraction with the dependence on column density removed 

using Equation (8) vs. T (bottom left panel). All data points used to make these plots passed the selection criteria described in 4.4. The 
color of each pixel is proportional to the logarithm of the number of data points located within the pixel. The solid lines show fits to the 
data (Section 6). The best-fit equations are listed on each plot in addition to the coefficient of determination (R^). 


of Figure 14 shows p vs. S. There is a clear anticorrela¬ 
tion between p and S {as = —0.67 ± 0.02, the coefficient 
of determination =0.47). We see no dependence of 
^ on iV (Figure 14, lower panel). 

We showed in Section 6 that the dependence of p on 
N can be removed using Equation (8) to create 
Similarly we can normalize out the dependence of p on 
S by calculating 


of with S is better than the correlation of p with 
S {R? =0.50 compared to 0.47). This indicates that 
both N and S contribute independently to the struc¬ 
ture seen in p. The fitted power-law exponents tend to 
be systematically shallower for the decorrelated and 
than for trends with p (see the first row of Table 
1 ), which might imply a weak underlying correlation be¬ 
tween N and S (see Figure 14, bottom panel). 


.[■5] 



(9) 


The top panel of Figure 15 shows that by removing 
the dependence of S from p the degree of correlation 
of with N increases {R? =0.35 compared to 0.30). 
Similarly, the bottom panel shows that the correlation 


7.2. Power-Law Fit p{N, S) 

As noted in Section 7.1, Figure 13 shows a color map 
of the median p binned two-dimensionally in S and N for 
ISRF-heated sightlines. The clear decrease of p with both 
increasing N and S is suggestive of a joint power-law rela¬ 
tionship. Here we derive a function p{N, S) that accounts 
for most of the structure seen in the p map. Specifically, 
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Table 1 

Power-law exponents of p vs. N and S 


Diffuse Emission 
Subtraction Method 

OLN 

as 

“p[Sl N 

"p[«i s 

Intermediate 

-0.58 ±0.02 

-0.67 ±0.02 

-0.46 ±0.01 

-0.58±0.01 

Conservative 

-0.53 ±0.02 

-0.67 ±0.02 

-0.39 ±0.01 

-0.56 ±0.01 

Aggressive 

-0.66 ±0.02 

-0.66 ±0.02 

-0.57±0.01 

-0.59±0.01 


Note. — The power-law exponents listed in this table are derived from linear fits 
of logp to log and log 5 as described in Sections 6 and 7.1. 
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Figure 13. Median p color-coded in bins of S and N for all 
ISRF-heated sightlines. The use of logarithmic scales for p, N, 
and S brings out the systematic relationship suggestive of power- 
laws. 


Table 2 

Fit parameters of p(N,S) from Equation (10) 


Diffuse Emission 
Subtraction Method 

aN 

as 

K 

Intermediate 

-0.45 ±0.01 

-0.60 ±0.01 

8.42 ±0.3 

Conservative 

-0.41 ±0.01 

-0.59 ±0.01 

6.92 ±0.3 

Aggressive 

-0.58 ±0.01 

-0.60 ±0.01 

10.98 ±0.3 


Note. — The power-law exponents (oat and as) and fitted 
constant (K) listed in this table are calculated from a two-variable 
power-law fit to N and S as described in Section 7.2. 


we adopt the joint power-law form 

\ogp{N,S) = K + aj\r log N + aslogS*, (10) 

where K, ajq and as are the free parameters. 

The exponents derived via a fit to Equation (10) are 
an =—0.45±0.01 and as =—0.60±0.01. Just as in 
Sections 6 and 7.1, errors in fit parameters are derived 
via bootstrapping (Table 2). We note that, as expected, 
aN and as derived from the two-variable power-law fit 
to N and S are identical within the error bars to ap[s]jv 
(the power-law fit to as a function of N) and ap[N]s 

(the power-law fit to as a function of 5), which were 

derived in Section 7.1 (also see Table 1). 

We can remove the dependence of p on V and S via 


[N,S] ^ PiP 


( 11 ) 
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Q. 
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Figure 14. Two-dimensional histograms showing correlations be¬ 
tween p, N, and S for ISRF-heated sightlines: p vs. S (top panel); 
and S vs N (bottom panel). All data points used to make these 
plots passed the selection criteria described in 4.4. The color is 
proportional to the logarithm of the number of data points located 
within each bin. The solid lines show fits to the data (Section 7.1). 

where is the decorrelated p for the fth data point, 

Pi is the measured polarization fraction for the fth data 
point, p{Ni,Si) is the value of the two-variable power- 
law fit for the fth data point, and p = 0.029 is the me¬ 
dian value of p. A comparison of the spatial distribution 
of p (middle panel) with (bottom panel) is shown 

in Figure 11. We discuss potential causes of residual 
structure in the plVS] Section 8.4. 

Finally, we quantify the degree to which our two- 
variable power-law fit p (V, S) can reproduce the ob¬ 
served dispersion in p. Figure 16 shows histograms of: 
our calculated logp (left panel); \ogp{N,S), the two- 
variable power-law fit values calculated for our N and 
S data points (center panel); and \ogp^^’^\ p with the 
derived dependence on N and S removed using Equa¬ 
tion (11) (right panel). For each of the three cases the 
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Figure 15. Two-dimensional histograms showing decorrelated 
p as a function of N and S for ISRF-heated sightlines: , the 

polarization fraction (p) with the dependence on S removed vs. 
N (top panel); and p^^^, p with the dependence on N removed vs. 
S (bottom panel). The color is proportional to the logarithm of 
the number of data points within the bin. The solid lines show the 
linear fits to the data (Section 7.1). 

median p is 0.029. Histograms of logj> rather than p are 
shown, because the fits are made in log space. The vari- 
ances of logp, logp {N, 5), and logpl^’'®! are 0.068, 0.045, 
and 0.023, respectively. Our model logp {N, S) repro¬ 
duces 66% of the variance in the logp map, which shows 
that our two-variable power-law fit model captures most 
of the physical effects that determine variations in frac¬ 
tional polarization in VelaC. 


7.3. Uncertainties in the Power-Law Fit Exponents 

The uncertainties of the exponents for the two-variable 
power-law fits ajq and as were estimated using the boot¬ 
strapping methods described in Section 6. However, as 
discussed in Section 3.5, the limiting uncertainty is our 
inability to precisely separate the contribution of back¬ 
ground/foreground dust from the polarized emission of 
VelaC. We repeated our analysis for maps of p and 
S made with the “conservative” and “aggressive” diffuse 
emission subtraction methods, to gauge the systematic 
uncertainty of our derived power-law exponents. Table 
1 gives power-law exponents derived from the individual 
correlations (Section 7.1) for the three different diffuse 
dust emission subtraction methods. Table 2 lists the ex¬ 
ponents ajv and as of the two-variable power-law fits 
again for all three diffuse emission subtraction methods. 
Systematic uncertainties relating to the choice of sub¬ 
traction method are seen to be ^0.1 for aN and ^0.01 


for as- 


8. DISCUSSION 

8 .1. Implieations of the Dependence of p on N and T 

In Section 6 we examined the dependence of the polar¬ 
ization fraction p on column density N and dust temper¬ 
ature T in VelaC. We divided our VelaC sightlines into 
two groups: those that show evidence of heating from 
the compact H II region RCW 36, and those sightlines 
where the temperature decreases as For the 

latter sightlines we suggested that the dust temperature 
is primarily set by exposure to the interstellar radiation 
field (ISRF), with high N sightlines having on average 
more shielding and therefore receiving less heating per 
unit mass from the ISRF. 

For the ISRF-heated sightlines, we find that p de¬ 
creases with increasing N and also that p increases as 
T increases. Depolarization for higher column density 
sightlines has been seen in many studies (see Section 
7.1). Vaillancourt & Matthews (2012) used the ratio 
of F(850/rm)/F(350^m) as a proxy for dust tempera¬ 
ture in two massive star forming clouds. They found 
that the polarization tended to decrease with increasing 
F'(850/rm)/F(350/rm), implying that warmer dust grain 
populations tend to have a higher p. This agrees with 
our result, but Vaillancourt & Matthews (2012) caution 
that variations in F(850^m)/F(350/rm) could be due to 
changes in dust spectral index, rather than just dust tem¬ 
perature. Our study is the first to fit p measurements 
within a molecular cloud as a function of both N and T. 
For the ISRF-heated sightlines, which are the majority 
of the sightlines, we find that the dependence of p on 
N is not separate from the dependence of p on T, since 
N and T are highly correlated. 

There are two general classes of explanations for our 
observations of p vs. N and T for the ISRF-headed sight¬ 
lines. We may have greater magnetic field disorder along 
high N (and therefore lower T) dust columns, or we may 
have a decrease in the intrinsic polarization efficiency 
(see Section 1) for such sightlines. In the first expla¬ 
nation the increased field disorder could arise because 
of a higher field disorder at high particle densities n, or 
because high N sightlines pass through more cloud ma¬ 
terial and therefore may sample different field directions 
at different locations along the line of sight (Jones 1989). 
Regarding the second possibile explanation for our ob¬ 
served p vs. N and T trends, note that in the radiative 
torques (RATs) model of grain alignment, “alignment 
torques” from an anisotropic radiation field are respon¬ 
sible for aligning the dust grain spin-axes with the local 
magnetic field (Lazarian & Hoang 2007). Grains near the 
surfaces of molecular clouds (low TV, high T) thus might 
be expected to show a higher average polarization frac¬ 
tion (Cho & Lazarian 2005). Alternatively, dust grain 
properties could change at high densities, e.g., grains 
could become rounder due to accretion of icy mantles 
(Whittet et al. 2008). 

For our RCW 36 heated sightlines there is a significant 
anticorrelation between p and N {R^ = 0.45). However, 
for these heated sightlines there is no correlation between 
N and T and no strong correlation between p and T 
{R^ = 0.03). This could indicate that the primary de¬ 
pendence of p is on A^, rather than T and that the cor- 
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Figure 16. Histograms of the logarithms of p before decorrelation (left panel), p{N,S) evaluated for N and S data using the model 
described in Equation (10) (center panel), and , the residual structure in p after normalizing out p (N, S) (right panel). The variance 

(cr^) of the distribution is given at the top of each panel. The quantity pA S] normalized so that the median p remained at 0.029 (see 
Equation (11)). 


relation of p with T only appears when there is a strong 
correlation of T with N. However, we caution that we 
have relatively few sightlines near RCW 36 (143 Nyquist- 
sampled sightlines compared to 2235 ISRF-heated sight¬ 
lines) so the lack of correlation between p and T could 
be caused by the angle of the magnetic field changing 
with respect to the line of sight, which would cause 
more spread in p. Indeed Figure 5 shows that near 
RCW 36 there are significant changes in the inferred mag¬ 
netic field orientation projected onto the plane of the sky. 

8 . 2 . Implications of the Two-variable Model p{N, S) 

In Section 7.2 we fit a model that describes p as a 
function with a power-law dependence on two variables, 
hydrogen column density N, and S, the dispersion in the 
polarization-angle on scales corresponding to our beam 
FWHM (2(5, or 0.5pc). The derived power-law expo¬ 
nents are aN = —0.45 ± 0.10 for the dependence on N, 
and as = —0.60 ± 0.01 for the dependence on S (see Sec¬ 
tion 7.3). Our p{N,S) fit is able to reproduce most of 
the structure seen in our logp maps. 

The decrease in p with increasing S can be attributed 
to changes in the magnetic field direction within the vol¬ 
ume probed by the BLASTPol beam. The mean mag¬ 
netic field orientation, $, is an average over both the 
beam area (0.5 pc) and along the length of the cloud in 
the line of sight direction, and is weighted by the density 
and intrinsic polarization efficiency (Section 1). Large 
values of S indicate a substantial change in $ on the 
scale of a beam, which implies a significant change in the 
orientation of polarization within the beam. This could 
be due to a sharp change in the magnetic field direction 
at some location within the cloud. Alternatively, it could 
indicate the overlap of two clouds, well separated along 
the line of sight, each with a different $. In either case we 
should see an overall decrease in the measured polariza¬ 
tion fraction, since some of the polarization components 
cancel. Planck Collaboration Int. XX (2015) note a de¬ 
crease of p with increasing S', both in their data and in 
corresponding MHD simulations (see Figure 19 of Planck 
Collaboration Int. XX 2015). However the Planck study 
sampled 5 x 10^° cm“^ < N < 10^^ cm“^ while our 
VelaC observations predominantly sample 10^^cm“^ < 
N < 10^^cm“^. Also direct comparison with their de¬ 
rived power-law exponents is difficult, since they fit S vs. 


p, thus minimizing the scatter in S, while we fit p vs. 
S, which minimizes scatter in p. Nevertheless they do 
find a significant anti-correlation of p and S in their data 
that is reproduced in their MHD simulations. In these 
simulations there is by contrast only a weak correlation 
of S with N (F. Levrier, private communication), just as 
we found in our data (Figure 14, lower panel). 

In Section 8.1, we discussed two classes of explana¬ 
tions for the observed p vs. N trend. The first class 
involves magnetic field disorder. An example is the work 
of Falceta-Gongalves et al. (2008). These authors were 
able to reproduce a decrease in p with increasing N via 
synthetic polarization maps made from supersonic, sub- 
Alfvenic MHD molecular cloud simulations, assuming 
uniform intrinsic polarization efficiency. Their power- 
law exponents ranged from 0 to —0.5, with models 
where the mean magnetic field was in the plane of the sky 
(7 = 0 °) having the steepest slope and models where the 
mean field was parallel to the line of sight (7 = 90°) hav¬ 
ing no dependence of p on N. In this theoretical study, 
the decrease in polarization for higher column density re¬ 
gions is due to an increase in the dispersion of the mag¬ 
netic field direction for high density regions. An analytic 
model by Jones (1989), is similarly able to reproduce a 
falling p vs. N for a medium having uniform intrinsic 
polarization efficiency. 

Our analysis shows only a weak correlation (or perhaps 
no correlation) between S and N (see Section 7). Thus 
on 0.5 pc scales, we find no significant increase in the dis¬ 
persion of $ for sightlines of increasing column density. 
Such an increase might be expected if disorder in the 
magnetic field direction increased in high density regions 
(for example due to accretion-driven turbulence as in 
Hennebelle & Andre 2013), or if the magnetic field were 
affected by large-scale gas motions near self gravitating 
filaments. In the above-mentioned theoretical study by 
Falceta-Gongalves et al. (2008), the authors showed that 
rarefied cloud regions show little variation in polarization 
direction whereas significant fluctuations in direction do 
occur within dense condensations. In this case, one might 
expect a positive correlation between S and N, which we 
do not see in our observations. 

The second class of explanations for the observed de¬ 
crease in p with increasing N involves changes in intrinsic 
polarization efficiency. This idea derives support from 
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the observations of Whittet et al. (2008). These authors 
measured the near-IR polarization of background stars 
in four nearby molecular clouds. For studies of polariza¬ 
tion of background starlight the quantity that is analo¬ 
gous to fractional polarization of dust emission is referred 
to as the “polarization efficiency”, defined as the frac¬ 
tional polarization of the starlight divided by the extinc¬ 
tion optical depth at the same wavelength px/t\. Whit¬ 
tet et al. (2008) found that the polarization efficiency 
in their clouds was consistent with a power-law depen¬ 
dence, V\It\ (X Because the inferred magnetic 

field direction is mostly uniform across the region stud¬ 
ied, they attributed all of the decrease in polarization 
efficiency with increasing N to changes in the intrinsic 
polarization efficiency. It is interesting to note that our 
power-law exponent (aAr = —0.45 ± 0.10) is similar to 
that found by Whittet et al. (2008). Other starlight po¬ 
larization studies have found power-law exponent values 
ranging from —0.34 to —1.0 (Goodman et al. 1995; Ger- 
akines et al. 1995; Arce et al. 1998; Ghapman et al. 2011; 
Alves et al. 2014; Gashman & Clemens 2014; Jones et al. 
2015). Ground-based studies of polarized thermal dust 
emission yield similar results. For example Matthews 
et al. (2002) examined p vs. Is 5 o for three clouds in Orion 
B South and found p oc (/gso/im)'^, with a ranging from 
-0.58 to -0.95. 

Which of the two general classes of explanations for 
the observed p vs. N trend best explains our observations 
of Vela C? Naively, the absence of a correlation between 
S and N would suggest that magnetic field disorder does 
not increase towards high N sightlines, which would im¬ 
ply that variation in intrinsic polarization efficiency is 
the more likely explanation. However, if the increased 
disorder in the field occurs on much smaller scales than 
0.5 pc, the scale probed by S, then S is not sensitive to 
the random component of the field and so we would not 
expect a correlation between N and S. We emphasize 
that detailed statistical comparisons with simulations of 
magnetized clouds that include variations in intrinsic po¬ 
larization efficiency are needed to fully understand the 
origin of the p vs. N anticorrelation (e.g., Soler et al. 
2013). Such comparisons are beyond the scope of the 
present paper. 

8.3. Analytic Models of p vs. N 

In Section 8.2 we advanced various explanations for 
the anticorrelation between p and N in VelaC. Here we 
consider an extreme case where all of the dependence of 
p on is due to reduced intrinsic polarization efficiency 
in shielded regions. Our goal is to quantify the ability 
of our measurements to trace magnetic fields deep inside 
the VelaC cloud under this pessimistic assumption. If 
most of the polarized emission comes from the outer dif¬ 
fuse layers of the cloud then our derived magnetic field 
orientations will not be sensitive to changes in the mag¬ 
netic field direction within dense structures embedded 
deep in the cloud. 

We model the efficiency of the dust along a given cloud 
sightline in emitting polarized radiation with e, where ep¬ 
silon is normalized such that e = ^ cos^ ( 7 ) IjAy, where 
f is the intrinsic polarization intensity as defined in Sec¬ 
tion 1 , Ay is the total dust extinction in the V band for 
that sightline, and 7 is the angle of the magnetic field 
with respect to the plane of the sky. For these mod- 


Line of Sight (z) 



Figure 17. Cartoon showing the parameterized depth into the 
cloud X for ^ slab model of a molecular cloud. The slab lies parallel 
to the plane of the sky. For a given position in the cloud along the 
line of sight (z) x is equal to the integrated visual dust extinction to 
the nearest cloud surface. The maximum value of x for a sightline 
of total visual extinction Ay is Ay!2. 
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Figure 18. BLASTPol measurements of the polarization frac¬ 
tion with the dependence of S normalized out vs. Ay. The 

solid line shows the results of a least-squares fit to the power-law 
decay intrinsic polarization efficiency (e) model with derived pa¬ 
rameters po = 0.038, Avcrit = 4.5mag and p = —1.21. The 
dashed-dotted line shows a fit to the constant e model with best-fit 
parameter po = 0.030. The dashed line shows a fit to the “skin 
depth” e model, where the best-fit parameters are po =0.039 
^Vcrit = 8.5 mag. 


els e = e(x), where x is the parameterized depth into 
the cloud, which is equal to the integrated visual extinc¬ 
tion to the nearest cloud surface as indicated in Figure 
17. Note that these models make a number of assump¬ 
tions: (a) that the cloud is isothermal; (b) that the dust 
emissivity does not change within the cloud; (c) that the 
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magnetic field direction is uniform; and (d) that the ge¬ 
ometry of the cloud is slab-like. 

As the p vs. N and p vs. S trends appear to be in¬ 
dependent (see Section 7) we compare predictions from 
our models with the polarization fraction decor- 

related from S (Figure 15, upper panel). Figure 18 
shows the predicted p from three models of e (y) com¬ 
pared to our measurements of vs. Ay, where here 
N has been converted to Ay assuming N{H)/Ay = 
2 X 10^^ cm“^ (Bohlin et al. 1978 with Ry =3.1). Be¬ 
cause of the normalization of pi"®! the overall level of po¬ 
larization is somewhat arbitrary, but of the right order. 

Here we describe the three plausible models of 
e (x) shown in Figure 18: 

Constant e Model: If the intrinsic polarization efficiency 
were constant throughout the cloud we would expect no 
dependence of p on N. The best fit to this model is 
shown as a dashed-dotted line in Figure 18. 

Skin Depth Model: Alternatively, we can consider a 
model where the intrinsic polarization efficiency is con¬ 
stant up to an extinction depth Xcrit and zero thereafter; 
i.e., a diffuse layer near the cloud surface is responsible 
for all of the polarized emission and the dust at cloud 
depths above Xcrit does not contribute to the polarized 
emission. For a sightline of total extinction Ay the max¬ 
imum value of X is Ay 12. We express e as 


e(x) = 


f eo , for X < Xcrit; 
I 0, for X > Xcrit, 


( 12 ) 



Figure 19. Models for the fraction of the total polarized intensity 
for a sightline of Ay = 40 mag (xmax = 20 mag), contributed 
by all dust at cloud depths < x' (see Equation (18)). The line 
color represents the power-law slope rj assumed: —0.8 (blue); and 
— 1.2 (black). Linestyle represents the Aycrit assumed: 3.0 mag 
(dashed); and 4.5 mag (solid). The red dotted line shows the ex¬ 
pected f(x') for dust of constant intrinsic polarization efficiency e. 


fractional polarization is then 



(i+i[»‘-i]) 


for a < 1; 
for a > 1. 


(17) 


where Eq is a constant, and from this we can calculate 
the total polarized intensity: 


fA-v/2 


P = 2 


' hi) dx 


f Ayeo, for Ay < Aycrit; 
(AvcritCo, for Ay > Aycrit, 


(13) 


where we define Aycrit = 2xcrit- The percentage polar¬ 
ization for a given sightline is then 

_ P _/P0, for Ay < Aycrit; . , 

^ ~ I “ \P0 , for Ay > Aycrit- ^ ’ 


In the “skin depth” model, p is constant for sightlines 
with Ay < Aycrit and decreases with a power-law slope 
of —1.0 for sightlines with Ay > Aycrit- The dashed line 
in Figure 18 shows a fit to the skin depth model. 
Power-law Model: Finally we consider a model where the 
polarization efficiency is constant up to Xcrit and there¬ 
after decreases as a power-law with coefficient t]: 


eix) 


eo, ^ for X < Xcrit; 

^0 (^) > for X > Xcrit- 


(15) 


This model simulates a constant e for the diffuse outer 
cloud layers and a decreasing e at greater cloud depths. 
The polarized intensity for a given sightline described by 
the power-law model is: 



for a < 1; 
, for a > 1, 


(16) 


where C = + 1 and a = Ay /Aycrit - The corresponding 


The power-law e model best fit parameters are po = 
0.038, Aycrit =4.5 mag and r] = —1.21 (Figure 18 solid 
line). Our power-law model fit would imply that at 
cloud depths of about two magnitudes or greater of visual 
extinction the intrinsic polarization efficiency decreases 
with depth into the cloud as ~ X~^'^^- 
It can be seen that both the skin-depth and power-law 
model capture the negative slope of the p^'^l vs. N curve 
for high N. For the purposes of quantifying our ability to 
trace magnetic fields deep within the cloud, we will use 
the power-law model as it seems to more closely follow 
the data points in Figure 18. We also caution that these 
are all simple models, so the fits should be taken merely 
as indicative of the trends of e with x- 
Using Equation (16) for a sightline of total dust ex¬ 
tinction Ay we can calculate the fractional contribution 
to the polarized intensity from cloud material at depths 
of < x' 


fix') 


Pjx!) 

-P(Xmax) 


(18) 


where by definition Xmax = Ay/2. Figure 19 shows 
/ (xO for a sightline of total dust extinction Ay = 
40 mag (about the largest found in Vela C) over the range 
of x' = 1 to x' = Xmax = 20 mag. Dashed and solid 
lines represent different assumptions for Aycrit, and line 
colors represent different power-law slopes g in Equa¬ 
tion (16). The solid black line is derived using the best- 
fit parameters. For comparison we also show the ex¬ 
pected behavior for the constant e model (red dotted 
line). For our best fit parameters 27% of the polar¬ 
ized emission comes from the outer 2.2 magnitudes of 
extinction, or the outer 2.2/20 =11% of the cloud. A 
further 47% of the total polarized emission comes from 
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Figure 20. Fraction of the dust column that is hidden, i.e. not 
traced by polarized emission as a function of Ay as described by 
Equation (21). The line color indicates rj and the linestyle indicates 
^Vcrit- 


2.2 mag < x < 10 mag, which accounts for 39% of the 
dust column. The most deeply embedded regions of the 
cloud (10 mag < x < 20 mag) contribute 50% of the 
total dust column but only 27% of the total polarized 
emission. Figure 19 also shows that steeper power-law 
slopes and lower Avcrit values would imply that more of 
the total polarized intensity measured comes from the 
outer diffuse cloud layers. 

To estimate the fraction of the cloud that, from the 
perspective of contributing to polarization, is “hidden” 
we first calculate the e weighted mean cloud depth (x): 


(x) 

Xedx 

C'" “‘X 

^Vcrit h (x) 

4 S (x) ’ 

(19) 

where g (x) 

— -^/(^O^Vcrit) 

(see Equation (16)) 

and 

h (x) is given 

by 



hix) = 

f 

for a < 1; 

(20) 

-b (a^ + 

^ — l) , for a > 1. 


If e were constant throughout the sightline then 
(x) would equal half of the maximum value of x giving 
(x) = 7lv/4. The fraction of the cloud that is hidden 
can then be roughly estimated as 


/hidden 1-0 


(x) 

Av/4’ 


( 21 ) 


which is shown in Figure 20. For %v = 10 mag (as¬ 
suming our best fit parameters) only 16% of the cloud is 
hidden. For a sightline of Ay = 40 mag about 48% of 
the cloud is hidden. 

In summary, for “moderate” dust column sightlines 
(yly < 10 mag) our polarization measurements sample 
most of the cloud (/hidden < 16%). So for sightlines 
with dust columns of Ay ~ 10 mag or less, the BLAST- 
Pol 500 /rm measurement of the magnetic field orienta¬ 
tion should be representative of the density-weighted av¬ 
erage magnetic field orientation along the sightline. For 
higher dust column sightlines, the fraction of the cloud 
that is not well sampled by our polarization measure¬ 
ments increases (/hidden 34% for Ay = 20 mag) and 


for our highest column sightlines (Ay ~ 40 mag) about 
half of the dust contributes little to the polarization mea¬ 
sured by BLASTPol. For these latter sightlines BLAST- 
Pol would not be sensitive to changes in magnetic field 
direction in the most deeply embedded cloud material. 
Recall that our model assumes that all of the decrease 
in p with N is due to lower intrinsic polarization effi¬ 
ciency of material deep within the cloud. If some of the 
decrease in p with N is due to increased field disorder 
along high N sightlines then the e drop-off with x would 
be shallower, which would decrease the fraction of the 
cloud that is hidden. 

As noted earlier, our model has many implicit assump¬ 
tions (slab-geometry, uniform dust temperature, power- 
law dependence of e(x))- In particular, the assumption 
of isothermality is clearly incorrect. Figure 10 shows that 
for the ISRF-heated sightlines included in this analysis 
the average temperature decreases with increasing col¬ 
umn density (and thus increasing Ay). For the temper¬ 
ature extremes of IIK and 15 K of the ISRF sightlines 
in Figure 10, we calculate that for the colder highest col¬ 
umn density sightlines the dust on average emits half as 
much radiation per unit mass at 500/rm compared with 
dust on the warmer more diffuse cloud sightlines. It is 
quite likely that the more deeply embedded dust grains 
in VelaC are colder, which implies they will contribute 
less than warmer grains near the cloud surfaces to both 
the total intensity and the measured polarized intensity. 
We therefore expect that the average magnetic field ori¬ 
entation inferred from polarization data will be weighted 
more towards the orientation in the warmer regions of the 
cloud. This will increase /hidden: even assuming uniform 
intrinsic polarization efficiency if the outer half of the 
dust grains had T = 15 K and the inner half of the dust 
grains had T = IIK then we find that /hidden = 20 %. 

To some degree this problem can be reduced by mea¬ 
suring polarization at millimeter wavelengths where the 
intensity of thermal dust emission is less sensitive to tem¬ 
perature. For detailed statistical comparisons of submil¬ 
limeter polarization data with synthetic observations of 
molecular clouds derived from numerical simulations it 
will be important to not only model the effects of grain 
alignment in simulation postprocessing but also include 
a realistic cloud temperature structure. Due to the afore¬ 
mentioned uncertainties that are related to the assump¬ 
tions in our model, our values of /hidden should be taken 
only as crude estimates. 

Despite these uncertainties, we note that, our re¬ 
sults are consistent with the findings of Cho & Lazarian 
(2005), who showed that dust grains can be aligned ef¬ 
ficiently by radiative torques at cloud depths x of up 
to 10 magnitudes in visual extinction. Bethell et al. 
(2007) found that the exact depth to which grains are 
aligned depends on grain size and on the degree of 
anisotropy of the local radiation field. Our model is also 
consistent with recent observations by Alves et al. (2014) 
who argue that their observations of submillimeter polar¬ 
ization of a starless core suggest loss of grain alignment at 
column densities higher than Ay = 30 mag. If e changes 
appreciably with x in the cloud then this might be re¬ 
vealed in the frequency dependence of p. Thus studying 
p at higher frequencies, as can be done using BLAST¬ 
Pol data, might provide further constraints. 
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8.4. Possible Causes of the Residual Structure in 

In Section 7.2 we showed that we can account for most 
of the variations in p that we observe in VelaC with 
a simple two-variable power-law fit p{N,S). Here, we 
consider a number of factors besides N and S which 
could contribute to the variance in p. The dispersion 
in the logarithm of the decorrelated fractional polariza¬ 
tion log is 0.15, which corresponds to a variance 

in Qf x.O X 10 “"^. 

If the variance in p were entirely due to the measure¬ 
ment uncertainty, then we would expect the variance in 
p to be described by: 


= = ( 22 ) 

where is the variance for each individual pi and n is 
the total number of data points. This value for 
is much smaller than the measured variance in ^ 

Measurement uncertainties thus play a minor part in the 
observed variance of 

A more likely possibility is that the variance in 
p[N,s] gggj^ jjj Figure 16 and the bottom panel of Figure 
11 is the result of changes in the direction of the mag¬ 
netic field with respect to the line of sight. The observed 
polarization of a population of dust grains aligned with 
respect to a uniform magnetic field depends on 7 , the 
angle between the magnetic field direction and the plane 
of the sky: 

P = Pma^xCOs'^l, (23) 

where Pmax is the polarization that would be observed 
if the magnetic field were parallel to the plane of the 
sky (7 = 0°). Our inferred magnetic field maps (Figure 
5) clearly show several large scale changes in magnetic 
field direction $. Corresponding large scale changes in 
7 would add width to the logp distribution. In theory a 
detailed statistical comparison of S and qj^ differ¬ 

ent angular scales could be used to gain insight into the 
three-dimensional structure of the magnetic field. How¬ 
ever, such a treatment is beyond the scope of the present 
paper. 

9. SUMMARY 

In this work we present 500 pm maps of the Vela C gi¬ 
ant molecular cloud from the 2012 flight of BLASTPol. 
Our polarization maps were calculated from Stokes I, 
Q and U maps with background/foreground diffuse po¬ 
larized emission subtracted as described in Section 3.5. 
These maps were used to calculate the inferred magnetic 
field orientation $ projected onto the plane of the sky. 
Overall we see a change in the magnetic field orienta¬ 
tion across the cloud, from perpendicular to the main 
cloud elongation direction in the south, to nearly par¬ 
allel to the cloud elongation in the north. We also see 
regions of sharp changes in the magnetic field direction, 
as traced by S, the average angular dispersion on scales 
corresponding to our beam (215 or 0.5pc scales). 

As a first step in our analysis of the VelaC data we ex¬ 
amine the dependence of polarization fraction p as a func¬ 
tion of column density N, dust temperature T, and local 
angular dispersion S for sightlines in four of the five cloud 


regions defined in Hill et al. (2011). The goal of this work 
is to look for empirical trends that can be compared to 
numerical simulations of molecular clouds. These trends 
can be used to learn about the magnetic field properties 
and intrinsic polarization efficiency within the cloud. As 
part of our analysis we separate our sightlines into those 
that appear to be primarily heated by the interstellar ra¬ 
diation field (ISRF) and the minority of sightlines that 
show evidence of heating from the compact H II region 
RCW36. 

Our main findings are as follows: 

1. For the ISRF-heated sightlines we find that p is 
anticorrelated with N and correlated with T, i.e., 
the polarization fraction decreases with increasing 
column density, and increases with increasing dust 
temperature. However, N and T are also highly 
anticorrelated with one another and normalizing 
out the power-law dependence of p with N re¬ 
moves the correlation with T. In the absence of 
bright internal sources it is expected that the den¬ 
sity structure of the cloud largely determines the 
observed T ; therefore we choose N as our indepen¬ 
dent variable in the subsequent analysis. For the 
RCW 36 heated sightlines where there is no cor¬ 
relation between N and T, we see no correlation 
between p and T but there is still a strong anticor¬ 
relation between p and N. This suggests that for 
the RCW 36-heated sightlines the important vari¬ 
able controlling p is N. 

2. We derive a two-variable power-law empiri¬ 
cal model p oc for the ISRF- 

heated sightlines, where aAr = —0.45 ± 0.10 and 
0:5 =—0.60 ± 0.01. This model can reproduce 
^ 66 % of the variance in logp. The decrease in 
p with increasing S is probably the result of changes 
in the magnetic field direction within the volume 
of the cloud sampled by the beam. The decrease 
in p with N could be caused by increased disor¬ 
der in the magnetic field for high column density 
sightlines or changes in the intrinsic polarization ef¬ 
ficiency (e.g., the fraction of aligned grains, or grain 
axis ratio) for deeply embedded cloud material. 

3. We do not find a strong correlation between 
N and S. This suggests that the disorder in the 
magnetic field does not increase significantly with 
density, which would in turn imply that the expla¬ 
nation for the decrease of p with increasing N is 
reduced intrinsic polarization efficiency for high 
N sightlines. However, this might not be the case. 
It might be that there is more disorder in the mag¬ 
netic field towards higher column density sightlines, 
but this disorder occurs on much smaller scales 
than 0.5 pc, the scale probed by S, such that S is 
not sensitive to the disordered magnetic field com¬ 
ponent. 

4. As a limiting case we consider the implications if 
the decrease in p with increasing N is due solely to 
reduced intrinsic polarization efficiency along high 
column density sightlines. In this case our BLAST¬ 
Pol measurements of the magnetic field orientation 
4> would preferentially sample the material closer 
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to the surface of the cloud and be less sensitive to 
changes in the field direction in the highly extin¬ 
guished regions deep within the cloud. We intro¬ 
duce a crude model in which the intrinsic polariza¬ 
tion efficiency is uniform in the outer cloud layers 
and then drops with a power-law dependence on 
the parameterized cloud depth y. From a fit of 
our observational data to this crude model we con¬ 
clude that for sightlines having Ay < 10 mag, $ is 
a reasonable measure of the average magnetic field 
direction along the line of sight, but for sightlines 
of Ay = 40 mag, much of the cloud (roughly the 
inner half) is not well traced by d*. This model 
might be a “worst-case” scenario because some of 
the decrease of p with Ay could arise from effects of 
magnetic field geometry not included in the model 
(e.g., more structure in the magnetic field along 
high column density sightlines). 

5. The remaining scatter in the polarization 

fraction with our derived power-law dependence on 
S and N normalized out, is too large to be ex¬ 
plained by our measurement uncertainties in p, but 
could be explained by variations in the angle of the 
magnetic field with respect to the plane of the sky. 

In this paper we have examined polarization trends for 
only one cloud. Other clouds with different properties, in 
particular different average angle 7 of the magnetic field 
with respect to the plane of the sky, might show differ¬ 
ent trends. To better constrain numerical simulations of 
star formation, our two-variable power-law fit should be 
repeated for a wide variety of clouds, which will presum¬ 
ably encompass a range of 7 values. Our study provides 
constraints for numerical simulations of molecular clouds; 
for at least one assumed value of 7 synthetic polarization 
observations of the simulations should be able to repro¬ 
duce (a) our two-variable power-law fit exponents and 
(b) the lack of correlation between N and S (on 0.5 pc 
scales). 
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APPENDIX 

A. DETAILED DESCRIPTION OF NULL TESTS 

We restrict consideration to the 250 pm maps for the 
purpose of the null tests because the larger number of 
detectors in this band allows coverage of the map area 
to remain complete even when splitting the data set in 
two. Furthermore, if the asymmetric beam shape is a 
significant source of systematics, these will be manifested 
more strongly at 250/rm, as the beam in this band is 
the least symmetric of the three. We would expect any 
regions that pass the null tests at 250/rm will also pass 
in the longer-wavelength bands. 

TODs were split into single raster scans, representing 
one complete raster of BLASTPol over the target map 
area at one half-wave plate position. Once the total data 
set was segmented using one of the four criteria described 
in Section 3.6, separate maps of Stokes /, Q, and U were 
made with TOAST (Section 3.4) for each of the two cate¬ 
gories. The diffuse emission removal described in Section 
3.5 was then performed for each map. For each null test 
criterion, we examined three metrics for evaluating sys¬ 
tematic disagreement between the data segments: 

1. Polarization fraction (p): Independent maps of p 
were produced using the /, Q, and U maps from 
each half of the data {pA and pb, where A and 
B generically represent the left and right sets, the 
early and late sets, etc.) A p residual map, Ap was 
calculated where Ap = {pA — pb)I the absolute 
value of which is absolute separation of each of 
PA and Pb from the mean p, {pA +pb)/‘^- The 
quantity Ap was taken to represent the uncertainty 
in p due to systematic sources of error, and we 
looked for regions in the full-data map where the 
calculated p is greater than 3Ap for each of the 4 
null tests described in Section 3.6. 

2. Polarization angle ('0): Analogously, two indepen¬ 
dent maps of 0 were calculated for each of the null 
tests (again, ip a and ipB, generically). Because 
a polarization measurement with 3a confidence in 
Q and U has an uncertainty in ip of about 10°, we 
looked for regions in the full-data map where the 
absolute difference between ip a and ipB was less 
than 20°. This standard is equivalent to requiring 
that the polarization-angle from each of the two 
data halves be consistent with the mean of ip a and 

i’B- 

3. Polarized intensity (P = + U"^): To examine 

systematic errors in P, we reproduced the proce¬ 
dure described in Section 3.4 of Matthews et al. 
(2014). Briefly, residual maps of Q and U were cal¬ 
culated as the difference between that parameter 
and its average value in the null test data halves. 
The Q and U residuals were then used to form a 
Pres = i/Qres + ^les- in the p metric described 
above. Pres was taken as the systematic uncertainty 
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in P, and we required the full-data measurement of 
P to be greater than 3Pres- 

B. POLARIZATION CONVENTIONS 

In this paper we discuss the polarized component of 
the dust emission (P) and the fractional polarization of 
dust emission (p), both of which can be derived from the 
linear polarization Stokes parameters: 

p = v/g2 + c/2, (Bi) 

and 

P=j- (B2) 

The associated angle of the polarization ip is 

■0 = ^ arctan {U, Q ), (B3) 

where the two argument arctan function computes 
arctan(g/P) while avoiding the ambiguity when Q = 0 
MJysr“^. The polarization angle 0 is defined from - 
90° to 90°. Our conventions for Q and U are such 
that 0° corresponds to North in Galactic coordinates 
and 0 increases East of Galactic North (counterclock¬ 
wise). This follows the lAU conventions (Hamaker & 
Bregman 1996), but differs from the HEALPix^^ conven¬ 
tion adopted for Planck data, where 0 increases West of 
Galactic North (clockwise). The apparent angle of the 
magnetic field projected on to the plane of the sky $ is 

$ = 0 -P (B4) 

It is important to note that $ is a tracer of the cloud 
magnetic field direction that is weighted by the efficiency 
of polarized dust emission averaged over the BLASTPol 
beam and along the line of sight. 

The variances of P, p and 0 are defined in Planck 
Gollaboration Int. XIX (2015) (Equations (B2)- (B4)). 
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